function [coef0, coef1, msglen, emin, cb, db] = embedimage(cmp_coef, img, t, k, qtable, msg, rands, typ, arg, varargin)
%EMBEDIMAGE embed message into JPEG image using MME.
%  [COEF0, COEF1, MSGLEN, EMIN, CB, DB] = EMBEDIMAGE(CMP_COEF, IMG, T, K, QTABLE, MSG, RANDS[, TYP, ARG])
%  The input argument IMG indicates the original JPEG coeffients. K is the
%  number of bits each segmentation has. MSG is the 0-1 stream to embed.
%  RANDS is the random stream to generate pusedo-random numbers. QTABLE is
%  an 8-by-8 matrix used in JPEG compression. T is the maximum number of
%  value to change.
%  CMP_COEF is just used for debugging, it is the coeffients read from JPEG
%  file.
%
%  The output argument COEF0 and COEF1 denotes the new JPEG DCT 
%  coeffients. LEN is the number of bytes of message embeded. 
%  EMIN is the overall minimal error got from IMG. CB is the bytes
%  of coeffients changed. DB is the bytes of coeffients in degeneration.
%  In original MMx (call GET_CHANGED_VALUE_ORIG), DB should be 0.
%
%  The randstream RANDS is used in the following order:
%  1. randperm       : rearrange the order of JPEG coeffients.
%  2. bitxor message : mess up with the message
%
%  Available TYP arguments are:
%    'bpc' : ARG specifies the maximum value of bpc, a value in [0..1].
%    'bpp' : ARG specifies the maximum value of bpp, a value in [0..1].
%    'bit' : ARG specifies the maximum number bits of message.
%
    global DEBUG;

%% parse input
    method = 1;
    typ = 'bit';
    arg = numel(msg);
    nvarargin = length(varargin);
    i = 1;
    while i <= nvarargin
        if strcmp(varargin{i}, 'Method')
            switch (varargin{i+1})
                case 'Original'
                    method = 0;
                case 'Modified'
                    method = 1;
                case 'Degenerate'
                    method = 2;
                case 'Adjust'
                    method = 3;
                otherwise
                    error('Unknown method.');
            end
        end
        i = i + 1;
    end
    %disp(method);
%% calculate errors
    n = 2^k - 1;
    %  O  -> c'           (non-round DCT coeffients)
    O = nonround_dct(img, qtable);
    % sym_o -> original symbol (original coeffients c'')
    sym_o = round(O);
    % r_o   -> R (c''-c') (rounding errors)
    % r_c -> S - C'       (rounding errors + embedding errors)
    % sym_c -> changed symbol (changed coeffients)
    % r_c = sym_c - O;
    switch (method)
        case 0
            r_o = sym_o - O;
            err = get_LSB_error(r_o, sym_o);
            sym_c = get_changed_value_orig(sym_o, r_o);            
        case 1
            r_o = sym_o - O;
            err = get_q_error(r_o, sym_o, qtable);
            sym_c = get_changed_value_orig(sym_o, r_o);
        case 2
            r_o = sym_o - O;
            err = get_LSB_error(r_o, sym_o);
            sym_c = get_changed_value(sym_o, r_o);
        case 3
            r_o = cmp_coef - O;
            err = get_LSB_error(r_o, cmp_coef);
            sym_c = get_changed_value_orig(cmp_coef, r_o);
        case 4
            r_o = sym_o - O;
            err = get_q_error2(r_o, sym_o, qtable);
            sym_c = get_changed_value_orig(sym_o, r_o);
    end
    clear img r_o r_c o;
    
%% DEBUG: check coeffients
    coef0 = sym_o;
    if DEBUG
        assert(numel(cmp_coef) == numel(coef0));    
        % The coeffients calculated from original image (COEF0) is different 
        % from the one from JPEG_READ (CMP_COEF)
        check_compressor(coef0, cmp_coef);
    end

%% skip DC
    [M0, N0] = size(sym_o);
    if method == 3
        sym_o = im2vec(cmp_coef,[8 8])';
    else
        sym_o = im2vec(sym_o,[8 8])';
    end
    sym_c = im2vec(sym_c,[8 8])';
    err = im2vec(err, [8 8])';
    ac_o = sym_o(:,2:end);
        
    ac_c = sym_c(:,2:end);
    ac_e = err(:,2:end);
   
%% random embed order
    [M, N] = size(ac_o);
    order = prand3(M*N, rands);    
    %order = prand1(M*N, rands); % DEBUG only, returned order is [1..M*N;1..M*N]
    ac_o = ac_o(order(1,:));
    ac_c = ac_c(order(1,:));
    ac_e = ac_e(order(1,:));
   
%% generate binary sequence from original message
    bin = arrayToBinary(msg);
    bin_len = numel(bin);
    if mod(bin_len, k) ~= 0
        % patch bin 
        old_len = bin_len;
        bin_len = ceil(bin_len / k) * k;
        bin = [bin zeros(1,bin_len - old_len)];
    end
    % bitxor
    bin = mod(bin(:) + round(rands.rand(bin_len,1)), 2);
    ind = find(ac_o ~= 0); % nonzero DCT
    ind_len = length(ind);

    % maximum embed rate
    if nargin < 9
        max_bit = numel(ac_o);
    else
        switch (typ)
            case 'bpc'
                if (arg < 0) || (arg > 1)
                    error('ARG should be a float value in [0..1].');
                end
                max_bit = ceil(ind_len * arg);
            case 'bpp'
                if (arg < 0) || (arg > 1)
                    error('ARG should be a float value in [0..1].');
                end
                max_bit = ceil(M0*N0*arg);
            case 'bit'
                if (arg < 0)
                    error('ARG should be a float value above 0.');
                end
                max_bit = ceil(arg);
            otherwise
                error('Unknown TYP value.');
        end
    end

    pc = ac_o;

%% embeding
    zlist = genzlist(t, k); % generate all possible Z out of loop to speedup
    etotal = 0; % total error
    cb = 0; % total change coeffients
    db = 0; % total degenerate coeffients
    msg_count = 1; % bits of message embeded
    i = 1; % number of used coeffients
    while (i <= ind_len)
        nexti = i + n;
        if (nexti - 1 > ind_len)
            % skip the last several coeffients & stop embedding ~
            w = sprintf('Not enough space to embed the whole message. %d/%d',msg_count,bin_len);
            warning(w);
            break;
        end
        
        %initial value
        current_msg = bin(msg_count:msg_count+k-1);
        current_coef = ac_o(ind(i:i+n-1));
        current_error = ac_e(ind(i:i+n-1));
        current_symbol = ac_c(ind(i:i+n-1));
        [nc, pos, emin] = mme_embedpart2(zlist, current_coef, current_error, current_symbol, current_msg, k);% embed
        degenerate_bits = (nc == 0);%ALWAYS TRUE: current_coef ~= 0
        degenerate = any(degenerate_bits);
        if DEBUG
            for j = i:i+n-1
                fprintf(1,'%d ',ind(j));
            end
            fprintf(1,'\t');
            debug_output(n, current_coef, nc, k, current_msg, degenerate, emin);
        end
        
        while degenerate
            % some Values change from 1 or -1 to 0
            % 1. set theses values to 0, change db
            num_degenerate = sum(degenerate_bits);

            if (nexti - 1 + num_degenerate > ind_len)
                w = sprintf('Not enough space to embed the whole message. %d/%d',msg_count,bin_len);
                warning(w);
                break;
            end
            db = db + num_degenerate;
            tmp = ac_o(ind(i:i+n-1));
            tmp(degenerate_bits) = 0;
            ac_o(ind(i:i+n-1)) = tmp;
            
            % 2. choose some more coefs (change:ind, nexti)
            new_ind = ind(i:i+n-1);
            new_ind = [new_ind(~degenerate_bits), ind(nexti:nexti+num_degenerate-1)];
            assert(numel(new_ind) == n);
            ind(i:i+n-1) = new_ind;
            nexti = nexti + num_degenerate;
            
            % 3. redo embed
            %current_msg = bin(msg_count:msg_count+k-1);
            current_coef = ac_o(ind(i:i+n-1));
            current_error = ac_e(ind(i:i+n-1));
            current_symbol = ac_c(ind(i:i+n-1));
            [nc, pos, emin] = mme_embedpart2(zlist, current_coef, current_error, current_symbol, current_msg, k);% embed
            degenerate_bits = (current_coef ~= 0) & (nc == 0);
            degenerate = any(degenerate_bits);
            if DEBUG
                for j = i:i+n-1
                    fprintf(1,'%d ',ind(j));
                end
                fprintf(1,'\t');
                debug_output(n, current_coef, nc, k, current_msg, degenerate, emin);
            end        
        end
        
        if degenerate
            % when only no more coeffients
            break;
        end
        
        ac_o(ind(i:i+n-1)) = nc;
        assert(isequal(pos, (nc ~= current_coef)));
        etotal = etotal + emin;
        msg_count = msg_count + k;
        cb = cb + sum(pos(:));
        clear pos emin;
        i = nexti;
        if (msg_count >= bin_len) || (msg_count >= max_bit) 
            % already embeded whole message or reach maximum rate
            break;
        end
    end
   
    ac_o = ac_o(order(2,:));% back to original order
    sym_o(:,2:end) = reshape(ac_o,[M N]);
    emin = etotal;
    coef1 = vec2im(sym_o',[0 0],[8 8],M0/8,N0/8);
    msglen = msg_count;

    assert(numel(coef0) == numel(coef1));
    %check_compressor(coef0, coef1);

    fprintf(1, 'image size = %d*%d px, total nonzero AC coef = %d\n', M0, N0, ind_len)
    fprintf(1, '%d / %d / %d of nonzero AC coef changed / degenerate / used \n',cb, db, i);
    fprintf(1, 'embeded message len = %db (%0.3f bpac, %0.3f bpp)\n', msg_count, msg_count / ind_len, msg_count / (M0*N0));
    fprintf(1, 'distortion error = %0.2f\n', emin);
end
